# Testes de hipótese #

# Artigo RAP: Políticos profissionais ou representantes de classe:
# Uma análise longitudinal da representação legislativa 
# no Brasil na democracia

# DataFrame: BD_Sainz-Albulquerque-Codato-2025

library(readxl)

BD_RAP_17082025 <- read_excel("BD_Sainz-Albulquerque-Codato-2025.xlsx")

View(BD_RAP_17082025)

# Carregar bibliotecas básicas #
library(dplyr)
library(tidyverse)
library(ggplot2)

# Estrutura do banco de dados

glimpse(BD_RAP_17082025)


######################################################
# Hipótese 1
######################################################

library(gmodels)
library(DescTools)
library(vcd)

# Filtrar base
dados_hip2 <- BD_RAP_17082025 |>
  filter(ocupacao_5_recod != "Outros/Nao informada") |>
  mutate(
    ideologia_recod = factor(ideologia_recod),
    ocupacao_5_recod = factor(ocupacao_5_recod),
    ano_eleicao_recod = as.numeric(ano_eleicao_recod)
  )

# Tabela de contingência
tabela_ocup_ideol <- table(dados_hip2$ocupacao_5_recod, dados_hip2$ideologia_recod)



# Qui-quadrado
teste_chi2 <- chisq.test(tabela_ocup_ideol)
teste_chi2

# V de Cramér
v_cramer <- CramerV(tabela_ocup_ideol)
v_cramer

# Visualizar resíduos ajustados
residuos_ajustados <- round(teste_chi2$stdres, 2)
residuos_ajustados

library(tidyr)

# Converter matriz de resíduos em data frame
residuos_df <- as.data.frame(as.table(residuos_ajustados)) |>
  rename(
    ocupacao = Var1,
    ideologia = Var2,
    residual = Freq
  )

residuos_df <- residuos_df |>
  mutate(
    ideologia = factor(ideologia, levels = c("Esquerda", "Centro", "Direita"))
  )

# Gráfico de calor
Mapa_H1 <- ggplot(residuos_df, aes(x = ideologia, y = ocupacao, fill = residual)) +
  geom_tile(color = "white", linewidth = 0.3) +
  geom_text(aes(label = round(residual, 2)), size = 5, color = "black") +
  scale_fill_gradient2(
    low = "#869B92", mid = "white", high = "#625A5D",
    midpoint = 0, limits = c(min(residuos_df$residual), max(residuos_df$residual)),
    name = "Resíduo"
  ) +
  labs(
    title = NULL,
    subtitle = NULL,
    x = "Ideologia",
    y = "Ocupação"
  ) +
  theme_minimal(base_size = 14) +
  theme(
    panel.grid = element_blank(),
    axis.text.x = element_text(angle = 0, vjust = 0.5, size = 16),
    axis.text.y = element_text(size = 16),
    legend.position = "right"
  )

Mapa_H1

ggsave(
  filename = "H1_mapa.png",
  plot = Mapa_H1,  # ou substitua por um objeto do gráfico, ex: plot = meu_grafico
  width = 12, height = 8, units = "in", dpi = 300
)

# Modelo Multinomial - Alternativa #

library(nnet)

# Criar subconjunto filtrado (se ainda não fez)
dados_multi <- BD_RAP_17082025 |>
  filter(ocupacao_5_recod != "Outros/Nao informada") |>
  mutate(
    ocupacao_5_recod = factor(ocupacao_5_recod),
    ideologia_recod = factor(ideologia_recod),
    ano_eleicao_recod = as.numeric(ano_eleicao_recod)
  )

dados_multi$ocupacao_5_recod <- relevel(dados_multi$ocupacao_5_recod, ref = "Politicos profissionais")

# Modelo básico: ocupação ~ ideologia + ano
modelo_multi <- multinom(ocupacao_5_recod ~ ideologia_recod + ano_eleicao_recod, data = dados_multi)

summary(modelo_multi)

exp(coef(modelo_multi))


# Modelo com interação

modelo_multi_inter <- multinom(ocupacao_5_recod ~ ideologia_recod * ano_eleicao_recod, data = dados_multi)

summary(modelo_multi_inter)

exp(coef(modelo_multi_inter))


# Visualização Modelo Multi

# Carregar pacotes necessários
library(dplyr)
library(tidyr)
library(ggplot2)
library(scales)

# 1. Construir grade com todas combinações possíveis
ideologias <- unique(dados_multi$ideologia_recod)
anos <- sort(unique(dados_multi$ano_eleicao_recod))
ocupacoes <- levels(dados_multi$ocupacao_5_recod)

grade_nova <- expand.grid(
  ideologia_recod = ideologias,
  ano_eleicao_recod = anos,
  ocupacao_5_recod = ocupacoes
)

# 2. Prever probabilidades com modelo multinomial
matriz_predita <- predict(
  modelo_multi,
  newdata = grade_nova,
  type = "probs"
)

# 3. Combinar grade com probabilidades preditas
df_preds <- as.data.frame(matriz_predita)
df_prob <- cbind(grade_nova, df_preds)

# 4. Transformar em formato longo para visualização
df_prob <- df_prob |>
  pivot_longer(
    cols = colnames(df_preds),
    names_to = "ocupacao_prevista",
    values_to = "probabilidade"
  )

# Ordenação

df_prob <- df_prob |>
  mutate(
    ideologia_recod = factor(ideologia_recod, levels = c("Esquerda", "Centro", "Direita"))
  )

# anos exatos das eleições
election_years <- sort(unique(dados_multi$ano_eleicao_recod))

Plot_H1 <- ggplot(df_prob,
       aes(x = as.numeric(ano_eleicao_recod),
           y = probabilidade,
           color = ideologia_recod,
           group = interaction(ideologia_recod, ocupacao_prevista))) +
  geom_line(linewidth = 0.5) +
  # se quiser, marque os pontos
  # geom_point(size = 1.6) +
  facet_wrap(~ ocupacao_prevista, scales = "free_y") +
  labs(
    title = NULL,
    x = "Ano da eleição", y = "Probabilidade ajustada", color = "Ideologia"
  ) +
  scale_y_continuous(labels = scales::percent_format(accuracy = 1)) +
  scale_x_continuous(breaks = election_years, labels = election_years,
                     minor_breaks = NULL,
                     expand = expansion(mult = c(0.01, 0.01))) +
  scale_color_manual(values = c("Esquerda"="#8E3B46","Centro"="#78C194","Direita"="#4C86A8")) +
  theme_minimal(base_size = 13) +
  theme(panel.grid.minor = element_blank(),
        panel.grid.major.x = element_blank(),
        panel.grid.major.y = element_blank(),
        axis.text.x = element_text(vjust = 0.5))

Plot_H1

ggsave(
  filename = "H1.png",
  plot = Plot_H1,  # ou substitua por um objeto do gráfico, ex: plot = meu_grafico
  width = 12, height = 8, units = "in", dpi = 300
)


# Extrair coeficientes e erros padrão
coefs <- summary(modelo_multi)$coefficients
std_errors <- summary(modelo_multi)$standard.errors

# Calcular z e p
z_values <- coefs / std_errors
p_values <- 2 * (1 - pnorm(abs(z_values)))

sig_stars <- function(p) {
  ifelse(p < 0.001, "***",
         ifelse(p < 0.01, "**",
                ifelse(p < 0.05, "*",
                       ifelse(p < 0.1, ".", ""))))
}

ocupacoes <- rownames(coefs)
variaveis <- colnames(coefs)

tabela_longa <- data.frame()

for (i in 1:nrow(coefs)) {
  for (j in 1:ncol(coefs)) {
    linha <- data.frame(
      Ocupacao = ocupacoes[i],
      Variavel = variaveis[j],
      Coeficiente = round(coefs[i, j], 3),
      Erro_Padrao = round(std_errors[i, j], 3),
      Valor_z = round(z_values[i, j], 2),
      p_valor = p_values[i, j],
      stringsAsFactors = FALSE
    )
    tabela_longa <- rbind(tabela_longa, linha)
  }
}

tabela_longa <- tabela_longa |>
  mutate(
    p_formatted = ifelse(p_valor < 0.001, "<0.001", 
                         ifelse(p_valor < 0.01, sprintf("%.3f", p_valor),
                                sprintf("%.3f", p_valor))),
    Significancia = sig_stars(p_valor),
    Coef_com_sig = paste0(sprintf("%.3f", Coeficiente), Significancia),
    EP_parenteses = paste0("(", sprintf("%.3f", Erro_Padrao), ")")
  ) |>
  mutate(
    Variavel_limpa = case_when(
      Variavel == "(Intercept)" ~ "Intercepto",
      Variavel == "ideologia_recodDireita" ~ "Direita",
      Variavel == "ideologia_recodEsquerda" ~ "Esquerda", 
      Variavel == "ano_eleicao_recod" ~ "Tempo (ano da eleição)",
      TRUE ~ Variavel
    ),
    Ocupacao_limpa = Ocupacao
  )

tabela_publicacao <- tabela_longa |>
  mutate(
    Ocupacao_display = ifelse(duplicated(Ocupacao_limpa), "", Ocupacao_limpa)
  ) |>
  select(Ocupacao_display, Variavel_limpa, Coef_com_sig, EP_parenteses, p_formatted) |>
  rename(
    "Ocupação" = Ocupacao_display,
    "Variável" = Variavel_limpa,
    "Coeficiente" = Coef_com_sig,
    "Erro Padrão" = EP_parenteses,
    "p-valor" = p_formatted
  )

print(tabela_publicacao)


tabela_compacta <- tabela_longa |>
  mutate(
    Resultado = paste0(sprintf("%.3f", Coeficiente), Significancia, 
                       " (", sprintf("%.3f", Erro_Padrao), ")"),
    Ocupacao_display = ifelse(duplicated(Ocupacao_limpa), "", Ocupacao_limpa)
  ) |>
  select(Ocupacao_display, Variavel_limpa, Resultado) |>
  rename(
    "Ocupação" = Ocupacao_display,
    "Variável" = Variavel_limpa,
    "Coef. (E.P.)" = Resultado
  )

print(tabela_compacta)

library(gt)

tabela_gt <- tabela_publicacao |>
  gt() |>
  tab_header(title = "Modelo Multinomial (sem interação)") |>
  cols_label(
    Ocupação = "Ocupação",
    Variável = "Variável",
    Coeficiente = "Coeficiente",
    `Erro Padrão` = "Erro padrão",
    `p-valor` = "Valor-p"
  ) |>
  cols_align(align = "left", columns = c(Ocupação, Variável)) |>
  cols_align(align = "center", columns = c(Coeficiente, `Erro Padrão`, `p-valor`)) |>
  opt_table_lines(extent = "none")

tabela_gt

# Hipótese 1 teste de efeitos marginais com robustez

library(marginaleffects)
library(dplyr)

## (1) AMEs para tempo
ame_tempo <- avg_slopes(
  modelo_multi_c,
  variables = "ano_c",
  type = "probs",
  cluster = ~ id_candidato
)

tabela_tempo <- as.data.frame(ame_tempo) |>
  mutate(
    ame_pp_ano   = 100 * estimate,
    ame_pp_ciclo = 4   * ame_pp_ano,
    IC95         = sprintf("[%.2f; %.2f]", 100*conf.low, 100*conf.high),
    `p-valor`    = signif(p.value, 3),
    Sig = ifelse(p.value < .001, "***",
                 ifelse(p.value < .01, "**",
                        ifelse(p.value < .05, "*",
                               ifelse(p.value < .1, ".", ""))))
  ) |>
  transmute(
    Ocupacao        = group,        # <- usar "group"
    `AME (%/ano)`   = round(ame_pp_ano, 2),
    `AME (%/ciclo)` = round(ame_pp_ciclo, 2),
    IC95, `p-valor`, Sig
  )

print(tabela_tempo)


library(marginaleffects)
library(dplyr)

# Comparações de ideologia (resumidas em médias por ocupação)
ame_ideologia <- avg_comparisons(
  modelo_multi_c,
  variables = "ideologia_recod",
  type = "probs"
)

tabela_ideologia <- as.data.frame(ame_ideologia) |>
  rename(Contraste = contrast) |>
  mutate(
    diff_pp  = 100 * estimate,
    IC95     = sprintf("[%.2f; %.2f]", 100*conf.low, 100*conf.high),
    `p-valor`= signif(p.value, 3),
    Sig      = ifelse(p.value < .001, "***",
                      ifelse(p.value < .01, "**",
                             ifelse(p.value < .05, "*",
                                    ifelse(p.value < .1, ".", ""))))
  ) |>
  transmute(
    Ocupacao = group,      # categoria de ocupação
    Contraste,             # comparação ideológica
    `Diferença (p.p.)` = round(diff_pp, 2),
    IC95, `p-valor`, Sig
  ) |>
  arrange(Ocupacao, Contraste)

print(tabela_ideologia)





######################################################
# Hipótese 2
######################################################

# Pacotes
library(dplyr)
library(tidyr)
library(ggplot2)
library(scales)
library(nnet)
library(DescTools)       
library(vcd)             
library(marginaleffects) 

# Modelo de regressão multinomial

dados_h3$ocupacao_5_recod <- relevel(dados_h3$ocupacao_5_recod, ref = "Profissionais liberais")

# sem interação
modelo_h3 <- multinom(ocupacao_5_recod ~ ideologia_recod + ano_c, data = dados_h3, trace = FALSE)

# com interação (conjunturas/tempo x ideologia)
modelo_h3_int <- multinom(ocupacao_5_recod ~ ideologia_recod * ano_c, data = dados_h3, trace = FALSE)

AIC(modelo_h3, modelo_h3_int)  # comparar ajuste
summary(modelo_h3_int)

# grid de previsão
anos <- sort(unique(dados_h3$ano_eleicao_recod))
grade_h3 <- expand.grid(
  ideologia_recod = levels(dados_h3$ideologia_recod),
  ano_eleicao_recod = anos
) |>
  mutate(ano_c = ano_eleicao_recod - 2010)

# prever (matriz com uma coluna por ocupação)
pred_mat <- predict(modelo_h3_int, newdata = grade_h3, type = "probs")
df_pred <- cbind(grade_h3, as.data.frame(pred_mat)) |>
  pivot_longer(
    cols = colnames(as.data.frame(pred_mat)),
    names_to = "ocupacao_prevista",
    values_to = "probabilidade"
  )

# 
df_pred <- df_pred |>
  mutate(
    ideologia_recod = factor(ideologia_recod, levels = c("Esquerda", "Centro", "Direita"))
  )

# gráfico
Plot_H2 <- ggplot(df_pred,
                  aes(x = ano_eleicao_recod, y = probabilidade,
                      color = ideologia_recod,
                      group = interaction(ideologia_recod, ocupacao_prevista))) +
  geom_line(linewidth = 0.5) +
  facet_wrap(~ ocupacao_prevista, scales = "free_y") +
  labs(
    title = NULL,
    x = "Ano da eleição", y = "Probabilidade ajustada", color = "Ideologia"
  ) +
  scale_y_continuous(labels = percent_format(accuracy = 1)) +
  scale_x_continuous(breaks = anos, labels = anos, minor_breaks = NULL) +
  scale_color_manual(values = c("Esquerda" = "#8E3B46", "Centro" = "#78C194", "Direita" = "#4C86A8")) +
  theme_minimal(base_size = 13) +
  theme(
    panel.grid.minor = element_blank(),
    panel.grid.major.y = element_blank(),
    panel.grid.major.x = element_blank(),
    axis.text.x = element_text(size = 12),
    axis.text.y = element_text(size = 12),
    axis.title.x = element_text(size = 12),
    strip.text = element_text(size = 14, face = NULL)  # AQUI aumenta o tamanho dos títulos
  )


Plot_H2

ggsave(
  filename = "H2.png",
  plot = Plot_H2,  # ou substitua por um objeto do gráfico, ex: plot = meu_grafico
  width = 12, height = 8, units = "in", dpi = 300
)

### ### ### ### ###

# Obter coeficientes e erros padrão
coefs <- summary(modelo_h3_int)$coefficients
ses <- summary(modelo_h3_int)$standard.errors

# Calcular valores z e p
z_vals <- coefs / ses
p_vals <- 2 * (1 - pnorm(abs(z_vals)))

# Combinar em um único data frame
tabela_resultados <- as.data.frame(
  cbind(
    coeficiente = as.vector(coefs),
    erro_padrao = as.vector(ses),
    z = as.vector(z_vals),
    p = as.vector(p_vals)
  )
)

# Adicionar colunas com nomes de variáveis e categorias
tabela_resultados$ocupacao <- rep(rownames(coefs), each = ncol(coefs))
tabela_resultados$variavel <- rep(colnames(coefs), times = nrow(coefs))

# Reorganizar colunas
tabela_resultados <- tabela_resultados |>
  select(ocupacao, variavel, coeficiente, erro_padrao, z, p)

# Exibir
print(tabela_resultados)

# Código para gerar tabela formatada do modelo multinomial H2

# Formatação da tabela do modelo multinomial para H2
# Baseado no código exato fornecido

library(dplyr)

# Função para estrelas de significância
sig_stars <- function(p) {
  ifelse(p < 0.001, "***",
         ifelse(p < 0.01, "**",
                ifelse(p < 0.05, "*",
                       ifelse(p < 0.1, ".", ""))))
}

# Extrair informações do modelo multinomial H2
summary_model <- summary(modelo_h3_int)
coefs <- summary_model$coefficients
std_errors <- summary_model$standard.errors

# Calcular valores z e p-valores
z_values <- coefs / std_errors
p_values <- 2 * (1 - pnorm(abs(z_values)))

# Criar tabela formatada
# As linhas representam as ocupações (Burocratas/Professores, Empresários, Trabalhadores)
# As colunas representam os preditores
# Categoria de referência: Profissionais liberais

ocupacoes <- rownames(coefs)
variaveis <- colnames(coefs)

# Criar data frame longo
tabela_longa <- data.frame()

for(i in 1:nrow(coefs)) {
  for(j in 1:ncol(coefs)) {
    linha <- data.frame(
      Ocupacao = ocupacoes[i],
      Variavel = variaveis[j],
      Coeficiente = round(coefs[i,j], 3),
      Erro_Padrao = round(std_errors[i,j], 3),
      Valor_z = round(z_values[i,j], 2),
      p_valor = p_values[i,j],
      stringsAsFactors = FALSE
    )
    tabela_longa <- rbind(tabela_longa, linha)
  }
}

# Formatar p-valores e adicionar estrelas
tabela_longa <- tabela_longa %>%
  mutate(
    p_formatted = ifelse(p_valor < 0.001, "<0.001", 
                         ifelse(p_valor < 0.01, sprintf("%.3f", p_valor),
                                sprintf("%.3f", p_valor))),
    Significancia = sig_stars(p_valor),
    Coef_com_sig = paste0(sprintf("%.3f", Coeficiente), Significancia),
    EP_parenteses = paste0("(", sprintf("%.3f", Erro_Padrao), ")")
  )

# Renomear variáveis para melhor apresentação
tabela_longa <- tabela_longa %>%
  mutate(
    Variavel_limpa = case_when(
      Variavel == "(Intercept)" ~ "Intercepto",
      Variavel == "ideologia_recodDireita" ~ "Direita",
      Variavel == "ideologia_recodEsquerda" ~ "Esquerda", 
      Variavel == "ano_c" ~ "Tempo (centrado em 2010)",
      Variavel == "ideologia_recodDireita:ano_c" ~ "Direita × Tempo",
      Variavel == "ideologia_recodEsquerda:ano_c" ~ "Esquerda × Tempo",
      TRUE ~ Variavel
    ),
    Ocupacao_limpa = case_when(
      Ocupacao == "Burocratas/Professores" ~ "Burocratas/Professores",
      Ocupacao == "Empresarios" ~ "Empresários",
      Ocupacao == "Trabalhadores" ~ "Trabalhadores",
      TRUE ~ Ocupacao
    )
  )

# Versão final para publicação
tabela_publicacao <- tabela_longa %>%
  mutate(
    Ocupacao_display = ifelse(duplicated(Ocupacao_limpa), "", Ocupacao_limpa)
  ) %>%
  select(Ocupacao_display, Variavel_limpa, Coef_com_sig, EP_parenteses, p_formatted) %>%
  rename(
    "Ocupação" = Ocupacao_display,
    "Variável" = Variavel_limpa,
    "Coeficiente" = Coef_com_sig,
    "Erro Padrão" = EP_parenteses,
    "p-valor" = p_formatted
  )

print(tabela_publicacao)

# Versão compacta (coeficiente e EP na mesma coluna)
cat("\n=== VERSÃO COMPACTA ===\n")
tabela_compacta <- tabela_longa %>%
  mutate(
    Resultado = paste0(sprintf("%.3f", Coeficiente), Significancia, 
                       " (", sprintf("%.3f", Erro_Padrao), ")"),
    Ocupacao_display = ifelse(duplicated(Ocupacao_limpa), "", Ocupacao_limpa)
  ) %>%
  select(Ocupacao_display, Variavel_limpa, Resultado) %>%
  rename(
    "Ocupação" = Ocupacao_display,
    "Variável" = Variavel_limpa,
    "Coef. (E.P.)" = Resultado
  )

print(tabela_compacta)

# Informações do modelo para nota de rodapé
cat("\n=== INFORMAÇÕES DO MODELO ===\n")
cat("Modelo: Regressão Logística Multinomial\n")
cat("Variável dependente: Ocupação (4 categorias)\n") 
cat("Categoria de referência: Profissionais liberais\n")
cat("Categoria de referência ideologia: Centro\n")
cat("Variável tempo: Centrada em 2010 (ano_c = ano_eleicao - 2010)\n")
cat(paste("N =", nrow(dados_h3), "observações\n"))
cat(paste("AIC modelo com interação =", round(AIC(modelo_h3_int), 1), "\n"))
cat(paste("AIC modelo sem interação =", round(AIC(modelo_h3), 1), "\n"))
cat("Significância: *** p<0.001; ** p<0.01; * p<0.05; . p<0.1\n")

library(gt)

# Criar a tabela formatada
tabela_gt <- tabela_publicacao |>
  gt() |>
  tab_header(title = "Modelo Multinomial com Interação") |>
  cols_label(
    Ocupação = "Ocupação",
    Variável = "Variável",
    Coeficiente = "Coeficiente",
    `Erro Padrão` = "Erro padrão",
    `p-valor` = "Valor-p"
  ) |>
  cols_align(align = "left", columns = c(Ocupação, Variável)) |>
  cols_align(align = "center", columns = c(Coeficiente, `Erro Padrão`, `p-valor`)) |>
  opt_table_lines(extent = "none")

# Visualizar no RStudio Viewer
tabela_gt

# Efeitos Marginais

library(dplyr)
library(marginaleffects)

# AME do tempo (ano_c) por cada ocupação (empilhado)
# categorias do desfecho (ocupações da H3)
cats <- levels(dados_h3$ocupacao_5_recod)

# função auxiliar que calcula AME para um outcome específico
ame_por_outcome <- function(out) {
  avg_slopes(
    modelo_h3_int,
    variables = "ano_c",
    type = "probs",
    outcome = out
  ) |>
    mutate(outcome = out)
}

# empilha resultados para todas as ocupações
ame_ano_h3 <- bind_rows(lapply(cats, ame_por_outcome))

# tabela limpa em pontos percentuais/ano e por ciclo (~4 anos)
ame_ano_tab <- ame_ano_h3 |>
  arrange(outcome) |>
  mutate(
    ame_pp_por_ano   = 100 * estimate,
    ame_pp_por_ciclo = 4   * ame_pp_por_ano
  ) |>
  select(outcome, ame_pp_por_ano, ame_pp_por_ciclo, conf.low, conf.high, p.value)

ame_ano_tab

# AME do tempo por ideologia (heterogeneidade), também por outcome

ame_por_outcome_by <- function(out) {
  avg_slopes(
    modelo_h3_int,
    variables = "ano_c",
    by = "ideologia_recod",
    type = "probs",
    outcome = out
  ) |>
    mutate(outcome = out)
}

ame_ano_by_ideol_h3 <- bind_rows(lapply(cats, ame_por_outcome_by))

ame_ano_by_tab <- ame_ano_by_ideol_h3 |>
  arrange(outcome, ideologia_recod) |>
  mutate(
    ame_pp_por_ano   = 100 * estimate,
    ame_pp_por_ciclo = 4   * ame_pp_por_ano
  ) |>
  select(outcome, ideologia_recod, ame_pp_por_ano, ame_pp_por_ciclo, conf.low, conf.high, p.value)

ame_ano_by_tab


# Todas as comparações pareadas (Direita–Centro, Esquerda–Centro, Esquerda–Direita)

comp_por_outcome <- function(out) {
  avg_comparisons(
    modelo_h3_int,
    variables = "ideologia_recod",
    comparison = "difference",
    type = "probs",
    outcome = out
  ) |>
    mutate(outcome = out)
}

ame_ideol_h3 <- bind_rows(lapply(cats, comp_por_outcome))

ame_ideol_tab <- ame_ideol_h3 |>
  arrange(outcome, contrast) |>
  mutate(
    diff_pp = 100 * estimate,
    li      = 100 * conf.low,
    ls      = 100 * conf.high
  ) |>
  select(outcome, contrast, diff_pp, li, ls, p.value)

ame_ideol_tab

# Por eleição
vc <- sandwich::vcovCL(modelo_h3_int, cluster = ~ ano_eleicao_recod)
avg_slopes(modelo_h3_int, variables = "ano_c", type = "probs", outcome = cats[1], vcov = vc)


# Tabelas

# 1) categorias do desfecho
cats <- levels(dados_h3$ocupacao_5_recod)

# 2) função auxiliar: AME de ano_c por outcome (coagir p/ data.frame aqui!)
ame_por_outcome <- function(out) {
  as.data.frame(
    avg_slopes(
      modelo_h3_int,
      variables = "ano_c",
      type = "probs",
      outcome = out
    )
  ) |>
    mutate(outcome = out)
}

# 3) empilhar resultados
ame_ano_h3 <- bind_rows(lapply(cats, ame_por_outcome))

# 4) formatar AME de tempo em p.p./ano e p.p./ciclo
sig_stars <- function(p) ifelse(p < .001,"***",
                                ifelse(p < .01,"**",
                                       ifelse(p < .05,"*",
                                              ifelse(p < .1,".",""))))

tabela_ame_tempo <- ame_ano_h3 |>
  # garantir nomes padronizados
  rename(p.value = any_of(c("p.value","Pr(>|z|)")),
         conf.low = any_of(c("conf.low","CI low")),
         conf.high = any_of(c("conf.high","CI high"))) |>
  mutate(
    ame_pp_por_ano   = 100 * estimate,
    ame_pp_por_ciclo = 4   * ame_pp_por_ano,
    IC95 = sprintf("[%.2f; %.2f]", 100*conf.low, 100*conf.high),
    stars = sig_stars(p.value)
  ) |>
  transmute(
    Outcome = outcome,
    `AME (%/ano)`   = round(ame_pp_por_ano, 2),
    `AME (%/ciclo)` = round(ame_pp_por_ciclo, 2),
    IC95,
    p.value = signif(p.value, 3),
    stars
  ) |>
  arrange(Outcome)

tabela_ame_tempo



# Tabela por ideologia

comp_por_outcome <- function(out) {
  as.data.frame(
    avg_comparisons(
      modelo_h3_int,
      variables  = "ideologia_recod",   # todas as comparações pareadas
      comparison = "difference",
      type = "probs",
      outcome = out
    )
  ) |>
    mutate(outcome = out)
}

ame_ideol_h3 <- bind_rows(lapply(cats, comp_por_outcome))

tabela_ame_ideol <- ame_ideol_h3 |>
  rename(p.value = any_of(c("p.value","Pr(>|z|)")),
         conf.low = any_of(c("conf.low","CI low")),
         conf.high = any_of(c("conf.high","CI high")),
         Contrast = any_of(c("contrast","Contrast"))) |>
  mutate(
    diff_pp = 100 * estimate,
    IC95 = sprintf("[%.2f; %.2f]", 100*conf.low, 100*conf.high),
    stars = sig_stars(p.value)
  ) |>
  transmute(
    Outcome = outcome,
    Comparação = Contrast,
    `Diferença (p.p.)` = round(diff_pp, 2),
    IC95,
    p.value = signif(p.value, 3),
    stars
  ) |>
  arrange(Outcome, Comparação)

tabela_ame_ideol



#### Tabela Clean ####

# estrelas de significância
sig_stars <- function(p) {
  ifelse(p < .001,"***",
         ifelse(p < .01,"**",
                ifelse(p < .05,"*",
                       ifelse(p < .1,".",""))))
}

# ordem amigável das ocupações (ajuste se quiser outra)
ord_ocup <- c("Empresarios", "Burocratas/Professores",
              "Profissionais liberais", "Trabalhadores")

# Tempo

tabela_ame_tempo_clean <-
  ame_ano_h3 |>
  # padroniza nomes que variam por versão
  rename(p.value  = any_of(c("p.value","Pr(>|z|)")),
         conf.low = any_of(c("conf.low","CI low")),
         conf.high= any_of(c("conf.high","CI high"))) |>
  # remove possíveis duplicatas exatas
  distinct(outcome, estimate, conf.low, conf.high, p.value, .keep_all = TRUE) |>
  mutate(
    ame_pp_ano   = 100 * estimate,
    ame_pp_ciclo = 4   * ame_pp_ano,
    IC95         = sprintf("[%.2f; %.2f]", 100*conf.low, 100*conf.high),
    stars        = sig_stars(p.value)
  ) |>
  transmute(
    Ocupacao        = outcome,
    `AME (%/ano)`   = round(ame_pp_ano, 2),
    `AME (%/ciclo)` = round(ame_pp_ciclo, 2),
    IC95,
    `p-valor`       = signif(p.value, 3),
    Sig             = stars
  ) |>
  arrange(match(Ocupacao, ord_ocup))

tabela_ame_tempo_clean


# Ideologia

tabela_ame_ideol_clean <-
  ame_ideol_h3 |>
  # padroniza nomes
  rename(p.value  = any_of(c("p.value","Pr(>|z|)")),
         conf.low = any_of(c("conf.low","CI low")),
         conf.high= any_of(c("conf.high","CI high")),
         Contrast = any_of(c("contrast","Contrast"))) |>
  # mantém só comparações de interesse
  filter(Contrast %in% c("mean(Direita) - mean(Centro)",
                         "mean(Esquerda) - mean(Centro)")) |>
  mutate(
    diff_pp = 100 * estimate,
    IC95    = sprintf("[%.2f; %.2f]", 100*conf.low, 100*conf.high),
    Sig     = sig_stars(p.value),
    Comparacao = case_when(
      Contrast == "mean(Direita) - mean(Centro)"  ~ "Direita − Centro",
      Contrast == "mean(Esquerda) - mean(Centro)" ~ "Esquerda − Centro",
      TRUE ~ Contrast
    )
  ) |>
  # remove duplicatas exatas
  distinct(outcome, Comparacao, diff_pp, IC95, p.value, .keep_all = TRUE) |>
  transmute(
    Ocupacao            = outcome,
    Comparacao,
    `Diferença (p.p.)`  = round(diff_pp, 2),
    IC95,
    `p-valor`           = signif(p.value, 3),
    Sig
  ) |>
  arrange(match(Ocupacao, ord_ocup),
          factor(Comparacao, levels = c("Direita − Centro","Esquerda − Centro")))

tabela_ame_ideol_clean


######################################################
# Hipótese 2 teste com erros padrao robustos cluster
######################################################

# === Pacotes ===
library(dplyr)
library(tidyr)
library(ggplot2)
library(scales)
library(nnet)
library(DescTools)       
library(vcd)             
library(marginaleffects) 
library(sandwich)
library(lmtest)
library(gt)

# ====================================================
# 1) Teste de associação
# ====================================================

# Base para H2: remover "Outros/Nao informada" e "Politicos profissionais"
dados_h3 <- BD_RAP_17082025 |>
  filter(!ocupacao_5_recod %in% c("Outros/Nao informada", "Politicos profissionais")) |>
  mutate(
    ocupacao_5_recod = droplevels(factor(ocupacao_5_recod)),
    ideologia_recod = factor(ideologia_recod, levels = c("Centro", "Direita", "Esquerda")),
    ano_eleicao_recod = as.numeric(ano_eleicao_recod),
    # centrar o ano para melhorar estabilidade numérica dos coeficientes
    ano_c = ano_eleicao_recod - 2010
  )

# Chi quadrado

tabela_h3 <- table(dados_h3$ocupacao_5_recod, dados_h3$ideologia_recod)

teste_chi2_h3 <- chisq.test(tabela_h3)
teste_chi2_h3
CramerV(tabela_h3)

residuos_h3 <- round(teste_chi2_h3$stdres, 2)
residuos_h3

# 1) Matriz de resíduos padronizados
res_mat <- round(teste_chi2_h3$stdres, 2)

# 2) Data frame “longo”
residuos_df <- as.data.frame(as.table(res_mat)) |>
  dplyr::rename(
    ocupacao  = Var1,
    ideologia = Var2,
    residual  = Freq
  ) |>
  dplyr::mutate(
    # garanta a ordem dos eixos (ajuste se quiser outra ordem)
    ideologia = factor(ideologia, levels = c("Centro", "Direita", "Esquerda")),
    ocupacao  = factor(ocupacao,  levels = c("Burocratas/Professores",
                                             "Empresarios",
                                             "Profissionais liberais",
                                             "Trabalhadores")),
    residual  = as.numeric(residual)
  )

# 3) Limites simétricos para a escala de cores
max_abs <- max(abs(residuos_df$residual), na.rm = TRUE)

# 4) Heatmap dos resíduos
ggplot(residuos_df, aes(x = ideologia, y = ocupacao, fill = residual)) +
  geom_tile(color = "white", linewidth = 0.3) +
  geom_text(aes(label = sprintf("%.2f", residual)), size = 4, color = "black") +
  scale_fill_gradient2(
    low = "#869B92", mid = "white", high = "#625A5D",
    midpoint = 0,
    limits = c(-max_abs, max_abs),
    name = "Resíduo\npadronizado"
  ) +
  labs(
    title = "Associação ocupação × ideologia (H3)",
    subtitle = "Resíduos padronizados ajustados do teste qui-quadrado",
    x = "Ideologia", y = "Ocupação"
  ) +
  theme_minimal(base_size = 14) +
  theme(
    panel.grid = element_blank(),
    axis.text.x = element_text(angle = 0, vjust = 0.5),
    legend.position = "right"
  )

# ====================================================
# 2) Modelos Multinomiais
# ====================================================
# Modelo de regressão multinomial

dados_h3$ocupacao_5_recod <- relevel(dados_h3$ocupacao_5_recod, ref = "Profissionais liberais")

# sem interação
modelo_h3 <- multinom(ocupacao_5_recod ~ ideologia_recod + ano_c, data = dados_h3, trace = FALSE)

# com interação (conjunturas/tempo x ideologia)
modelo_h3_int <- multinom(ocupacao_5_recod ~ ideologia_recod * ano_c, data = dados_h3, trace = FALSE)

AIC(modelo_h3, modelo_h3_int)  # comparar ajuste
summary(modelo_h3_int)

# grid de previsão
anos <- sort(unique(dados_h3$ano_eleicao_recod))
grade_h3 <- expand.grid(
  ideologia_recod = levels(dados_h3$ideologia_recod),
  ano_eleicao_recod = anos
) |>
  mutate(ano_c = ano_eleicao_recod - 2010)

# prever (matriz com uma coluna por ocupação)
pred_mat <- predict(modelo_h3_int, newdata = grade_h3, type = "probs")
df_pred <- cbind(grade_h3, as.data.frame(pred_mat)) |>
  pivot_longer(
    cols = colnames(as.data.frame(pred_mat)),
    names_to = "ocupacao_prevista",
    values_to = "probabilidade"
  )

# 
df_pred <- df_pred |>
  mutate(
    ideologia_recod = factor(ideologia_recod, levels = c("Esquerda", "Centro", "Direita"))
  )

# gráfico
Plot_H2 <- ggplot(df_pred,
                  aes(x = ano_eleicao_recod, y = probabilidade,
                      color = ideologia_recod,
                      group = interaction(ideologia_recod, ocupacao_prevista))) +
  geom_line(linewidth = 0.5) +
  facet_wrap(~ ocupacao_prevista, scales = "free_y") +
  labs(
    title = NULL,
    x = "Ano da eleição", y = "Probabilidade ajustada", color = "Ideologia"
  ) +
  scale_y_continuous(labels = percent_format(accuracy = 1)) +
  scale_x_continuous(breaks = anos, labels = anos, minor_breaks = NULL) +
  scale_color_manual(values = c("Esquerda" = "#8E3B46", "Centro" = "#78C194", "Direita" = "#4C86A8")) +
  theme_minimal(base_size = 13) +
  theme(
    panel.grid.minor = element_blank(),
    panel.grid.major.y = element_blank(),
    panel.grid.major.x = element_blank(),
    axis.text.x = element_text(size = 12),
    axis.text.y = element_text(size = 12),
    axis.title.x = element_text(size = 12),
    strip.text = element_text(size = 14, face = NULL)  # AQUI aumenta o tamanho dos títulos
  )


Plot_H2

ggsave(
  filename = "H2.png",
  plot = Plot_H2,  # ou substitua por um objeto do gráfico, ex: plot = meu_grafico
  width = 12, height = 8, units = "in", dpi = 300
)

# ====================================================
# 4) Probabilidades preditas (gráfico H2)
# ====================================================
anos <- sort(unique(dados_h2$ano_eleicao_recod))
grade_h2 <- expand.grid(
  ideologia_recod   = levels(dados_h2$ideologia_recod),
  ano_eleicao_recod = anos
) |>
  mutate(ano_c = ano_eleicao_recod - 2010)

pred_mat <- predict(mod_h2_int, newdata = grade_h2, type = "probs")
df_pred <- cbind(grade_h2, as.data.frame(pred_mat)) |>
  pivot_longer(cols = colnames(as.data.frame(pred_mat)),
               names_to = "ocupacao_prevista", values_to = "probabilidade")

Plot_H2 <- ggplot(df_pred, aes(x = ano_eleicao_recod, y = probabilidade,
                               color = ideologia_recod,
                               group = interaction(ideologia_recod, ocupacao_prevista))) +
  geom_line(linewidth = 0.6) +
  facet_wrap(~ ocupacao_prevista, scales = "free_y") +
  scale_y_continuous(labels = percent_format(accuracy = 1)) +
  scale_x_continuous(breaks = anos) +
  scale_color_manual(values = c("Esquerda" = "#8E3B46", "Centro" = "#78C194", "Direita" = "#4C86A8")) +
  labs(x = "Ano da eleição", y = "Probabilidade ajustada", color = "Ideologia") +
  theme_minimal(base_size = 13)

print(Plot_H2)

ggsave("H2.png", Plot_H2, width = 12, height = 8, dpi = 300)


# ====================================================
# Hipótese 2 - AME por ocupação com SE clusterizado
# ====================================================

library(dplyr)
library(marginaleffects)

# 1) Lista de categorias do desfecho
cats <- levels(dados_h2$ocupacao_5_recod)

# 2) Função auxiliar: calcula AME de ano_c por ocupação
ame_por_outcome <- function(out) {
  as.data.frame(
    avg_slopes(
      mod_h2_int,
      variables = "ano_c",
      type = "probs",
      outcome = out,
      cluster = ~ id_candidato   # <- correção: cluster por candidato
    )
  ) |>
    filter(term == "ano_c") |>    # mantém só o efeito marginal principal
    mutate(outcome = out)
}

# 3) Empilhar resultados
ame_ano_cluster <- bind_rows(lapply(cats, ame_por_outcome))

# 4) Função de estrelas de significância
sig_stars <- function(p) {
  ifelse(p < .001,"***",
         ifelse(p < .01,"**",
                ifelse(p < .05,"*",
                       ifelse(p < .1,".",""))))
}

# 5) Tabela formatada para publicação
tabela_ame_tempo_clean <- ame_ano_cluster |>
  rename(p.value  = any_of(c("p.value","Pr(>|z|)")),
         conf.low = any_of(c("conf.low","CI low")),
         conf.high= any_of(c("conf.high","CI high"))) |>
  mutate(
    ame_pp_ano   = 100 * estimate,
    ame_pp_ciclo = 4 * ame_pp_ano,
    IC95         = sprintf("[%.2f; %.2f]", 100*conf.low, 100*conf.high),
    Sig          = sig_stars(p.value)
  ) |>
  transmute(
    Ocupacao        = outcome,
    `AME (%/ano)`   = round(ame_pp_ano, 2),
    `AME (%/ciclo)` = round(ame_pp_ciclo, 2),
    IC95,
    `p-valor`       = signif(p.value, 3),
    Sig
  )

# 6) Visualizar resultado
print(tabela_ame_tempo_clean)

# ====================================================
# Hipótese 2 - AME por ocupação e ideologia (SE clusterizado)
# ====================================================

# 1) Lista de categorias do desfecho
cats <- levels(dados_h2$ocupacao_5_recod)

# 2) Função auxiliar: calcula AME de ano_c por ocupação × ideologia
ame_por_outcome_by <- function(out) {
  as.data.frame(
    avg_slopes(
      mod_h2_int,
      variables = "ano_c",
      type = "probs",
      outcome = out,
      by = "ideologia_recod",     # <- agora estratificado por ideologia
      cluster = ~ id_candidato    # <- correção: cluster por candidato
    )
  ) |>
    filter(term == "ano_c") |>    
    mutate(outcome = out)
}

# 3) Empilhar resultados
ame_ano_cluster_by <- bind_rows(lapply(cats, ame_por_outcome_by))

# 4) Função de estrelas de significância
sig_stars <- function(p) {
  ifelse(p < .001,"***",
         ifelse(p < .01,"**",
                ifelse(p < .05,"*",
                       ifelse(p < .1,".",""))))
}

# 5) Tabela formatada para publicação
tabela_ame_tempo_by <- ame_ano_cluster_by |>
  rename(p.value  = any_of(c("p.value","Pr(>|z|)")),
         conf.low = any_of(c("conf.low","CI low")),
         conf.high= any_of(c("conf.high","CI high")),
         Ideologia = ideologia_recod) |>
  mutate(
    ame_pp_ano   = 100 * estimate,
    ame_pp_ciclo = 4 * ame_pp_ano,
    IC95         = sprintf("[%.2f; %.2f]", 100*conf.low, 100*conf.high),
    Sig          = sig_stars(p.value)
  ) |>
  transmute(
    Ocupacao        = outcome,
    Ideologia,
    `AME (%/ano)`   = round(ame_pp_ano, 2),
    `AME (%/ciclo)` = round(ame_pp_ciclo, 2),
    IC95,
    `p-valor`       = signif(p.value, 3),
    Sig
  ) |>
  arrange(Ocupacao, Ideologia)

# 6) Visualizar resultado
print(tabela_ame_tempo_by)



######################################################
# Hipótese 3
######################################################

library(dplyr)
library(nnet)

# Base H4: sem "Outros/Nao informada" e sem "Politicos profissionais"
dados_h4 <- BD_RAP_17082025 |>
  filter(!ocupacao_5_recod %in% c("Outros/Nao informada", "Politicos profissionais")) |>
  mutate(
    ocupacao_5_recod = droplevels(factor(ocupacao_5_recod)),
    ideologia_recod  = factor(ideologia_recod, levels = c("Centro","Direita","Esquerda")),
    ano_eleicao_recod = as.numeric(ano_eleicao_recod),
    ano_c = ano_eleicao_recod - 2010  # centraliza para estabilidade numérica
  )


# Modelo

# Modelo reduzido (base H3): sem interação tempo×ideologia
mod_h3_reduz <- multinom(
  ocupacao_5_recod ~ ideologia_recod + ano_c,
  data  = dados_h4,
  trace = FALSE
)

# Modelo completo (H4): COM interação tempo×ideologia
mod_h4_compl <- multinom(
  ocupacao_5_recod ~ ideologia_recod * ano_c,
  data  = dados_h4,
  trace = FALSE
)

# --- LRT manual (robusto a classe 'multinom') ---
ll_red   <- logLik(mod_h3_reduz)
ll_full  <- logLik(mod_h4_compl)
lr_stat  <- 2 * (as.numeric(ll_full) - as.numeric(ll_red))
df_diff  <- attr(ll_full, "df") - attr(ll_red, "df")
p_lrt    <- pchisq(lr_stat, df = df_diff, lower.tail = FALSE)

cat("=== LRT Multinomial (H4) ===\n")
cat("LogLik reduzido:", round(as.numeric(ll_red), 2), "\n")
cat("LogLik completo:", round(as.numeric(ll_full), 2), "\n")
cat("LR stat:", round(lr_stat, 2), " | df:", df_diff, " | p-valor:", format.pval(p_lrt), "\n")


# Log Linear

# Tabela 3D
tab3d_h4 <- xtabs(~ ocupacao_5_recod + ideologia_recod + ano_eleicao_recod, data = dados_h4)

# Modelo com TODAS as interações de 2ª ordem (sem 3ª)
m2 <- loglin(tab3d_h4, margin = list(
  c("ocupacao_5_recod","ideologia_recod"),
  c("ocupacao_5_recod","ano_eleicao_recod"),
  c("ideologia_recod","ano_eleicao_recod")
), fit = TRUE)

# Modelo COMPLETO com interação de 3ª ordem
m3 <- loglin(tab3d_h4, margin = list(
  c("ocupacao_5_recod","ideologia_recod","ano_eleicao_recod")
), fit = TRUE)

# Diferença de deviances (LRT)
delta_dev <- m2$lrt - m3$lrt
delta_df  <- m2$df  - m3$df
p_ll      <- pchisq(delta_dev, df = delta_df, lower.tail = FALSE)

cat("\n=== Log-linear (H4) ===\n")
cat("ΔDeviance:", round(delta_dev, 2), " | Δdf:", delta_df, " | p-valor:", format.pval(p_ll), "\n")


# Visualização

# 0) Sair de debug se estiver no prompt Browse[1]>
# (faça manualmente no console: digite Q e tecle Enter)

# --- pressupõe que 'dados_h4' e 'mod_h4_compl' já existem como no seu script ---

library(dplyr)
library(tidyr)
library(ggplot2)
library(scales)
library(nnet)

# 1) Grid de previsão (tempo como numérico contínuo)
anos <- sort(unique(dados_h4$ano_eleicao_recod))
grid_pred <- expand.grid(
  ideologia_recod   = levels(dados_h4$ideologia_recod),
  ano_eleicao_recod = anos
) |>
  mutate(
    ideologia_recod   = factor(ideologia_recod, levels = c("Esquerda","Centro","Direita")),
    ano_eleicao_recod = as.numeric(ano_eleicao_recod),
    ano_c             = ano_eleicao_recod - 2010
  )

# 2) Probabilidades previstas por ocupação (multinom::predict)
#    Retorna matriz n x K (K = categorias de ocupacao)
prob_mat <- predict(mod_h4_compl, newdata = grid_pred, type = "probs")

# 3) Para longo: uma linha por (ideologia, ano, ocupacao)
df_pred <- cbind(grid_pred, as.data.frame(prob_mat)) |>
  pivot_longer(
    cols = setdiff(colnames(as.data.frame(prob_mat)), character(0)),
    names_to = "outcome",
    values_to = "estimate"
  ) |>
  mutate(
    outcome = factor(outcome, levels = levels(dados_h4$ocupacao_5_recod))
  ) |>
  arrange(outcome, ideologia_recod, ano_eleicao_recod)

# 4) Gráfico (sem ribbon, para evitar IC incorreto)
ggplot(
  df_pred,
  aes(x = ano_eleicao_recod, y = estimate,
      color = ideologia_recod, group = ideologia_recod)
) +
  geom_line(linewidth = 0.5) +
  facet_wrap(~ outcome, scales = "free_y") +
  scale_y_continuous(labels = percent_format(accuracy = 1)) +
  scale_x_continuous(
    breaks = anos, labels = anos, minor_breaks = NULL,
    expand = expansion(mult = c(0.01, 0.01))
  ) +
  scale_color_manual(values = c("Esquerda"="#8E3B46","Centro"="#78C194","Direita"="#4C86A8")) +
  labs(
    title = NULL,
    x = "Ano da eleição", y = "Probabilidade ajustada",
    color = "Ideologia"
  ) +
  theme_minimal(base_size = 13) +
  theme(
    panel.grid.minor = element_blank(),
    panel.grid.major.x = element_blank(),
    panel.grid.major.y = element_blank()
  )

# ====================================================
# Hipótese 3 - AMEs por ideologia com SE clusterizados
# ====================================================

# Função auxiliar: AME do tempo (ano_c) por ocupação × ideologia
ame_por_outcome_by <- function(out) {
  as.data.frame(
    avg_slopes(
      mod_h4_compl,
      variables = "ano_c",
      type = "probs",
      outcome = out,
      by = "ideologia_recod",   # estratificação por ideologia
      cluster = ~ id_candidato  # robustez clusterizada
    )
  ) |>
    filter(term == "ano_c") |>    
    mutate(outcome = out)
}

# Empilhar resultados
ame_h3_cluster_by <- bind_rows(lapply(cats, ame_por_outcome_by))

# Tabela formatada
tabela_ame_h3_by <- ame_h3_cluster_by |>
  rename(
    p.value   = any_of(c("p.value","Pr(>|z|)")),
    conf.low  = any_of(c("conf.low","CI low")),
    conf.high = any_of(c("conf.high","CI high")),
    Ideologia = ideologia_recod
  ) |>
  mutate(
    ame_pp_ano   = 100 * estimate,
    ame_pp_ciclo = 4 * ame_pp_ano,
    IC95         = sprintf("[%.2f; %.2f]", 100*conf.low, 100*conf.high),
    Sig          = ifelse(p.value < .001,"***",
                          ifelse(p.value < .01,"**",
                                 ifelse(p.value < .05,"*",
                                        ifelse(p.value < .1,".",""))))
  ) |>
  transmute(
    Ocupacao        = outcome,
    Ideologia,
    `AME (%/ano)`   = round(ame_pp_ano, 2),
    `AME (%/ciclo)` = round(ame_pp_ciclo, 2),
    IC95,
    `p-valor`       = signif(p.value, 3),
    Sig
  ) |>
  arrange(Ocupacao, Ideologia)

# Visualizar resultado
print(tabela_ame_h3_by)


######################################################
# Hipótese 4
######################################################

# Pacotes
library(dplyr)
library(tidyr)
library(broom)
library(DescTools)   # BreslowDayTest

# ----------------------------
# 1) Preparação dos dados
# ----------------------------
dados_h5 <- BD_RAP_17082025 |>
  filter(!ocupacao_5_recod %in% c("Outros/Nao informada", "Politicos profissionais")) |>
  mutate(
    ideologia_recod   = factor(ideologia_recod, levels = c("Centro","Direita","Esquerda")),
    ano_eleicao_recod = as.numeric(ano_eleicao_recod)
  )

# Escolha do ponto de corte do "tempo" (ajuste se quiser)
# Sugestão: 2010 separa pré/pós conjunturas centrais
cut_ano <- 2010

dados_h5 <- dados_h5 |>
  mutate(
    tempo_bin = ifelse(ano_eleicao_recod > cut_ano, "pos", "pre"),
    tempo_bin = factor(tempo_bin, levels = c("pre","pos"))
  )

# Ocupações-alvo (ajuste a lista conforme seu desenho)
ocupacoes_alvo <- c("Empresarios", "Burocratas/Professores",
                    "Profissionais liberais", "Trabalhadores")

# ----------------------------
# 2) Funções de teste por ocupação
# ----------------------------
monta_array_2x2xk <- function(df, ocupacao_nome) {
  df |>
    mutate(
      occ_bin = factor(ifelse(ocupacao_5_recod == ocupacao_nome, "alvo","outros"),
                       levels = c("outros","alvo"))
    ) |>
    xtabs(~ occ_bin + tempo_bin + ideologia_recod, data = _)
}

testa_h5_para <- function(ocupacao_nome, df = dados_h5) {
  arr <- monta_array_2x2xk(df, ocupacao_nome)
  
  # MH para 2x2xK: OR comum e teste de associação estratificada
  mh <- mantelhaen.test(arr)  # base R
  
  # Homogeneidade dos ORs entre estratos (ideologias)
  bd <- DescTools::BreslowDayTest(arr)  # requer 2x2xk
  
  tibble(
    ocupacao      = ocupacao_nome,
    or_mh         = as.numeric(mh$estimate),                      # OR comum (MH)
    ci_low        = as.numeric(mh$conf.int[1]),
    ci_high       = as.numeric(mh$conf.int[2]),
    p_mh          = mh$p.value,                                   # p do MH
    p_breslowday  = bd$p.value                                    # p da homogeneidade
  )
}

# ----------------------------
# 3) Rodar para todas as ocupações
# ----------------------------
resultados_h5 <- bind_rows(lapply(ocupacoes_alvo, testa_h5_para)) |>
  mutate(
    interpretacao = case_when(
      p_breslowday < 0.05 ~ "ORs heterogêneos entre blocos (diferenças entre esquerda/centro/direita)",
      TRUE                ~ "Homogeneidade não rejeitada (OR comum plausível)"
    )
  )

resultados_h5





#### Revisão Hipótese 1 do Artigo RAP ####

library(nnet)
library(dplyr)
library(tidyr)
library(ggplot2)
library(scales)
library(gt)

# Preparar base com ano centralizado
dados_multi <- BD_RAP_17082025 |>
  filter(ocupacao_5_recod != "Outros/Nao informada") |>
  mutate(
    ocupacao_5_recod = factor(ocupacao_5_recod),
    ideologia_recod = factor(ideologia_recod),
    ano_eleicao_recod = as.numeric(ano_eleicao_recod),
    ano_c = ano_eleicao_recod - 2010  # Centralizar o ano em 2010
  )

# Definir categoria de referência
dados_multi$ocupacao_5_recod <- relevel(dados_multi$ocupacao_5_recod, ref = "Politicos profissionais")

modelo_multi <- multinom(ocupacao_5_recod ~ ideologia_recod + ano_c, data = dados_multi)

# Coeficientes e erros padrão
coefs <- summary(modelo_multi)$coefficients
std_errors <- summary(modelo_multi)$standard.errors

# Estatísticas
z_values <- coefs / std_errors
p_values <- 2 * (1 - pnorm(abs(z_values)))

# Função para asteriscos de significância
sig_stars <- function(p) {
  ifelse(p < 0.001, "***",
         ifelse(p < 0.01, "**",
                ifelse(p < 0.05, "*",
                       ifelse(p < 0.1, ".", ""))))
}

# Montar tabela longa
ocupacoes <- rownames(coefs)
variaveis <- colnames(coefs)
tabela_longa <- data.frame()

for (i in 1:nrow(coefs)) {
  for (j in 1:ncol(coefs)) {
    tabela_longa <- rbind(tabela_longa, data.frame(
      Ocupacao = ocupacoes[i],
      Variavel = variaveis[j],
      Coeficiente = round(coefs[i, j], 3),
      Erro_Padrao = round(std_errors[i, j], 3),
      Valor_z = round(z_values[i, j], 2),
      p_valor = p_values[i, j],
      stringsAsFactors = FALSE
    ))
  }
}

tabela_longa <- tabela_longa |>
  mutate(
    p_formatted = ifelse(p_valor < 0.001, "<0.001", sprintf("%.3f", p_valor)),
    Significancia = sig_stars(p_valor),
    Coef_com_sig = paste0(sprintf("%.3f", Coeficiente), Significancia),
    EP_parenteses = paste0("(", sprintf("%.3f", Erro_Padrao), ")"),
    Variavel_limpa = case_when(
      Variavel == "(Intercept)" ~ "Intercepto",
      Variavel == "ideologia_recodDireita" ~ "Direita",
      Variavel == "ideologia_recodEsquerda" ~ "Esquerda",
      Variavel == "ano_c" ~ "Tempo (centrado em 2010)",
      TRUE ~ Variavel
    ),
    Ocupacao_limpa = Ocupacao
  )

# Tabela padrão
tabela_publicacao <- tabela_longa |>
  mutate(Ocupacao_display = ifelse(duplicated(Ocupacao_limpa), "", Ocupacao_limpa)) |>
  select(Ocupacao_display, Variavel_limpa, Coef_com_sig, EP_parenteses, p_formatted) |>
  rename(
    "Ocupação" = Ocupacao_display,
    "Variável" = Variavel_limpa,
    "Coeficiente" = Coef_com_sig,
    "Erro Padrão" = EP_parenteses,
    "p-valor" = p_formatted
  )

# Tabela compacta
tabela_compacta <- tabela_longa |>
  mutate(
    Resultado = paste0(sprintf("%.3f", Coeficiente), Significancia,
                       " ", EP_parenteses),
    Ocupacao_display = ifelse(duplicated(Ocupacao_limpa), "", Ocupacao_limpa)
  ) |>
  select(Ocupacao_display, Variavel_limpa, Resultado) |>
  rename(
    "Ocupação" = Ocupacao_display,
    "Variável" = Variavel_limpa,
    "Coef. (E.P.)" = Resultado
  )

# Visualizar tabela formatada
tabela_gt <- tabela_publicacao |>
  gt() |>
  tab_header(title = "Hipótese 1: Modelo de Regressão Logística Multinomial") |>
  cols_label(
    Ocupação = "Ocupação",
    Variável = "Variável",
    Coeficiente = "Coeficiente",
    `Erro Padrão` = "Erro padrão",
    `p-valor` = "Valor-p"
  ) |>
  cols_align(align = "left", columns = c(Ocupação, Variável)) |>
  cols_align(align = "center", columns = c(Coeficiente, `Erro Padrão`, `p-valor`)) |>
  opt_table_lines(extent = "none")

# Exibir
tabela_gt
